library(Glimma)
library(limma)
library(GlimmaV2)
data(lymphomaRNAseq)
rnaseq <- lymphomaRNAseq
# 7 samples
# each instance is a gene
# attribute is number of counts per sample
# data.frame(rnaseq$counts)
# add lane
groups <- data.frame(genotype=rnaseq$samples$group,
lane= as.character(c(rep(4,5),3,3)),
miscCont=c(rep(4000,5),300,250),
miscDisc=c("blue","red",rep("green",5)))
# add libsize
groups <- rnaseq$samples$group
# fit
design <- model.matrix(~0+groups)
contrasts <- cbind(Smchd1null.vs.WT=c(-1,1))
# convert raw counts to logCPM values by automatically extracting libsizes and normalisation factors from x
vm <- voomWithQualityWeights(rnaseq, design=design)
fit <- lmFit(vm, design=design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
dt <- decideTests(fit)
summary(dt)
## Smchd1null.vs.WT
## Down 1478
## NotSig 24424
## Up 1277
testing <- function(fit, counts, groups=1:ncol(counts), logTransform=FALSE)
{
if (is.null(colnames(counts)))
{
samples <- seq_len(ncol(counts))
}
else
{
samples <- colnames(counts)
}
# process counts
cols <- colnames(counts)
rownames(counts) <- make.names(cols)
if (logTransform)
{
counts <- as.matrix(edgeR::cpm(counts, log=TRUE))
}
else
{
counts <- as.matrix(counts)
}
counts <- t(counts)
return(counts)
if (!is.numeric(groups))
{
# reorder samples to group levels
groups <- factor(groups)
counts <- counts[order(groups), ]
sample.cols <- sample.cols[order(groups)]
samples <- samples[order(groups)]
groups <- sort(groups)
}
expData <- data.frame(
Sample = samples,
cols = as.hexcol(sample.cols),
Group = groups,
counts)
expData
}
# object has a bunch of genes
nrow(fit$genes)
## [1] 27179
# fold change is the attribute here
# each instance is a gene;
# the col vector names are the geneIDs
# the col vector attribute is fold change (Smchd1null vs WT)
# data.frame(fit$coefficients)
# log cpm is the attribute here
# averaged across all arrays in the original linear model fit
# each instance is a gene;
# the col vector names are the geneIDs
# the col vector attribute is mean expression
# data.frame(fit$Amean)
# for now assume fit object
# this is where info for the MD scatter plot comes from!
# data.frame(vm$E)
# fit$design
GlimmaV2("MA", fit)
## Warning in readLines(con): incomplete final line found on 'C:/Users/hkari/
## OneDrive/onenote/R/win-library/3.6/GlimmaV2/htmlwidgets/GlimmaV2.yaml'